home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / STIRLING.C < prev   
C/C++ Source or Header  |  1995-12-01  |  5KB  |  228 lines
  1. /* ============ */
  2. /* stirling.c    */
  3. /* ============ */
  4. #include <defcodes.h>
  5. #include <math.h>
  6. #include <mconf.h>
  7. /* ------------------- */
  8. /* FUNCTION PROTOTYPES */
  9. /* ------------------- */
  10. # undef F
  11. # if defined(__STDC__) || defined(__PROTO__)
  12. #    define    F( P )    P
  13. # else
  14. #    define    F( P )    ()
  15. # endif
  16. /* INDENT OFF */
  17. extern    double    BinCoef F((double, double));
  18. extern    LDBL    Stirling1 F((double, double));
  19. extern    LDBL    Stirling2 F((double, double));
  20.  
  21. # undef F
  22. /* INDENT ON */
  23. #define    MAX_ELEMS    100        /* Max Recursion List */
  24. #define MAXNUM    1.701411834604692317316873e38;    /* 2**127 */
  25.  
  26. /* ==================================================================== */
  27. /* Stirling1 - Returns Stirling Number of the First Kind for S(n,m)    */
  28. /* ==================================================================== */
  29. LDBL
  30. Stirling1(double n, double m)
  31. {
  32.     LDBL    RetVal;
  33.  
  34.     /* ---------------------------------------------------------------- */
  35.     /* The expressions used are from Handbook of Mathmetical Functions,    */
  36.     /* ed. by M. Abramowitz and I. A. Stegun, U.S. Government Printing    */
  37.     /* Office, 1964, Paragraph 24.1.3 and Table 24.3.            */
  38.     /* ---------------------------------------------------------------- */
  39.     /* -------------------------------------------------------- */
  40.     /*    S(n, m)   = S(n-1,m-1) + (n-1) * S(n-1, m), n >= m >= 1 */
  41.     /*    S(n, n)   = 1                        */
  42.     /*    S(n, n-1) = [n*(n-1)]/2                    */
  43.     /*    S(n, 0)   = 0                        */
  44.     /*    S(n, 1)   = (n-1)!                    */
  45.     /* -------------------------------------------------------- */
  46.  
  47.     if (m < 0 || n < 0 || m > n || n-m > MAX_ELEMS)
  48.     {
  49.     char    Cause[32] = "Stirling1(): ";
  50.  
  51.     strcat(Cause,
  52.         (m < 0) ? "(m < 0)" :
  53.         (n < 0) ? "(n < 0)" :
  54.         (m > n) ? "(m > n)" : "(n-m > 100)");
  55.  
  56.     mtherr(Cause, DOMAIN);
  57.  
  58.     RetVal = MAXNUM;
  59.     }
  60.     else if (m == n)
  61.     {
  62.     RetVal = 1;
  63.     }
  64.     else if (m == (n-1))
  65.     {
  66.     RetVal = (m*n)/2;
  67.     }
  68.     else if (m == 0)
  69.     {
  70.     RetVal = 0;
  71.     }
  72.     else if (m == 1)
  73.     {
  74.     RetVal = (n < 2) ? 1 : fac((int)n-1);
  75.     }
  76.     else
  77.     {
  78.     LDBL    List[MAX_ELEMS], *P, Val;
  79.     int    k;
  80.  
  81.     /* ------------------------------- */
  82.     /* Set up recursion run for m = 1  */
  83.     /* (these are factorials (k-1)! ). */
  84.     /* ------------------------------- */
  85.     P = List;
  86.     *P = 1;
  87.  
  88.     for (k = 2; k <= n - m; ++k)
  89.     {
  90.         Val = k * (*P++);
  91.  
  92.         *P = Val;
  93.     }
  94.     /* -------------------------- */
  95.     /* Recurse runs from 2 to m.  */
  96.     /* There are (n-m-1) elements */
  97.     /* per run.              */
  98.     /* -------------------------- */
  99.     for (k = 2; k <= m; ++k)
  100.     {
  101.        int     j;
  102.        LDBL    Mult;
  103.  
  104.         P = List;
  105.  
  106.         /* ------------------------- */
  107.         /* First element is k(k+1)/2 */
  108.         /* ------------------------- */
  109.         *P = (k * (k+1)) >> 1;
  110.         Mult = k + 1;
  111.  
  112.         for (j = 2; j <= n-m; ++j)
  113.         {
  114.         Val = (Mult++) * *P++;
  115.  
  116.         *P += Val;
  117.         }
  118.     }
  119.  
  120.     RetVal = *P;
  121.     }
  122.  
  123.     return (RetVal);
  124. }
  125. /* ==================================================================== */
  126. /* Stirling2 - Returns Stirling Number of the Second Kind for S(m,n)    */
  127. /* ==================================================================== */
  128. LDBL
  129. Stirling2(double n, double m)
  130. {
  131.     long double  RetVal;
  132.  
  133.     /* ---------------------------------------------------------------- */
  134.     /* The expressions used are from Handbook of Mathmetical Functions,    */
  135.     /* ed. by M. Abramowitz and I. A. Stegun, U.S. Government Printing    */
  136.     /* Office, 1964, Paragraph 24.1.4 and Table 24.4.            */
  137.     /* ---------------------------------------------------------------- */
  138.     /* -------------------------------------------------------- */
  139.     /*    S(n, m)   = S(n-1,m-1) + m * S(n-1, m), n >= m >= 1     */
  140.     /*    S(n, n)   = 1                        */
  141.     /*    S(n, n-1) = [n*(n-1)]/2                    */
  142.     /*    S(n, 0)   = 0                        */
  143.     /*    S(n, 1)   = 1                        */
  144.     /*    S(n, 2)   = 2^n - 1                    */
  145.     /* -------------------------------------------------------- */
  146.  
  147.     if (m < 0 || n < 0 || m > n || n-m > MAX_ELEMS)
  148.     {
  149.     char    Cause[32] = "Stirling2(): ";
  150.  
  151.     strcat(Cause,
  152.         (m < 0) ? "(m < 0)" :
  153.         (n < 0) ? "(n < 0)" :
  154.         (m > n) ? "(m > n)" : "(n-m > 100)");
  155.  
  156.     mtherr(Cause, DOMAIN);
  157.  
  158.     RetVal = MAXNUM;
  159.     }
  160.     else if (m == n)
  161.     {
  162.     RetVal = 1;
  163.     }
  164.     else if (m == (n-1))
  165.     {
  166.     RetVal = (m*n)/2;
  167.     }
  168.     else if (m == 0)
  169.     {
  170.     RetVal = 0;
  171.     }
  172.     else if (m == 1)
  173.     {
  174.     RetVal = 1;
  175.     }
  176.     else if (m == 2)
  177.     {
  178.     RetVal = ldexp(0.5, (int)(n)) - 1;
  179.     }
  180.     else
  181.     {
  182.     LDBL    List[MAX_ELEMS], *P, Val;
  183.     int    k;
  184.  
  185.     /* ------------------------------- */
  186.     /* Set up recursion run for m = 1  */
  187.     /* (these are Pk = 2^(k+1) - 1 ).  */
  188.     /* ------------------------------- */
  189.     P = List;
  190.  
  191.     for (k = 1; k <= n - m; ++k)
  192.     {
  193.         *P++ = ldexp(2., k) - 1;
  194.     }
  195.  
  196.     /* -------------------------- */
  197.     /* Recurse runs from 3 to m.  */
  198.     /* There are (n-m-1) elements */
  199.     /* per run.              */
  200.     /* -------------------------- */
  201.     for (k = 3; k <= m; ++k)
  202.     {
  203.        int     j;
  204.        LDBL    Mult;
  205.  
  206.         P = List;
  207.  
  208.         /* ------------------------- */
  209.         /* First element is k(k+1)/2 */
  210.         /* ------------------------- */
  211.         *P = (k * (k+1)) >> 1;
  212.  
  213.         Mult = k;
  214.  
  215.         for (j = 2; j <= n-m; ++j)
  216.         {
  217.         Val = Mult * *P++;
  218.  
  219.         *P += Val;
  220.         }
  221.     }
  222.  
  223.     RetVal = *P;
  224.     }
  225.  
  226.     return (RetVal);
  227. }
  228.